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r this program takes in a raw 8 bit image, with x and y dimensions as inputs also, 
and encodes that image with a specified number, spitting the encoded image out */ 

T the current version, started just after sign_it.sav3, now adds the three parameters 
of dc, root of the gray value, and a constant times the local gradient */ \J : . 



#include "pinecone.h" " 
#include "image.h" 

#define NOISE_MAX 10000000 

#define BITS 16 

#define HALFJVIAX 32768 . 

#define MAX 65536 

#define MAX1 65535 

#define GV.POWER 0.4 

#define GRADJTHRESHOLD 50.0 

#deflne DISPLAY JT 1 

unsigned char *img; 
long bit[BITS]; 
long gidjgv.gidjgrad; 
long xdim.ydim; 
float *gv,*grad; 
Colorindex *disp; 
Colorindex *scale_disp; 
Colorindex *pdisp; 

int helplines = 2; 
char *helpfl = 
{ 

"usage: signjt filename [xdim(default 510) ydim(default 256)]\n", 
"\nMust input both xdim and ydim if not using the defaults " 

); 

int fiil_bits(int code_number){ 
long i.tmp; 

printf("The binary value "); 

tmp = HALF_MAX; 
for(i=0;i<16;i++){ 

if(code_number>=tmp){ 
printffr'); 

code_number-=tmp; 
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bit[Q=1; 

} 

else { 

printf("0 M ); 
btt[i]=0; 

} 

tmp/=2; 

} 

printff will be wallpapered across the image\n\n"); 
return(O); 

} 

int load_sca1e_disp(void){ 
long i j; 

pdisp = scale_disp+20; 
for(p0y<15;j++){ 

for(i=1;i<20;i++){ 

*pdisp « (unsigned char)255; 

pdisp += 20; 

} 

pdisp += 20; 

} 

pdisp = scale_disp+20+(35*400); 
forG=0;j<15;j++){ 

for(i=1;i<20;i++){ 

*pdisp - (unsigned char)255; 
pdisp += 20; 

} 

pdisp += 20; 

} 

return(O); 
) 

int draw_scale(float value Jong which){ 
long i; 
float flemp; 

value *=5; 

memcpytdisp.scale^disp^OOOO); 
ftemp = 20.0 * value; 
if(ftemp > 399.0)ftemp=399.0; 
pdisp = &disp[1 5*400 + (int)ftemp]; 
for(i=0;i<20;i++){ 

*pdisp = (unsigned char)255; 



pdisp +- 400; 

} 

winset(which); 
rectwritefO.O.SQMS.disp); 

return(O); 
} 

int calculate jgvO{ 
unsigned char *pimg; 
long i; 
float *pgv; 

pimg = img; 
pgv = gv; 

for(i=0;i<(ydim*xdim);i++){ 

*(pgv++) B P owf( (float)*(pimg++) , GV POWER ); 

} 

retum(O); 
} 

int calculate_grad(){ 

unsigned char *pimg1 l *pimg2,*pimg; 

long 

float avg,*pgrad; 

pimg = &img[1+xdim]; 
pimgl = &img[1]; 
pimg2 = &img[1 +2*xdim]; 

pgrad = grad; 
for(i=0;i<(xdim*ydim);i++)*(pgrad++) = 0.0; 

pgrad = &grad[1+xdim]; 
for(i=1 ;i<(ydim-1);i++){ 

for0=1 j<(xdim-1);j++){ 
avg = (floatr(pimg1-1)+ (float)*(pimg1++)+ (fioat)*pimg1 + (float)*(pimg- 
1)+(float)*(pimg+1)+(float)*(pimg2-1)+(float)*(pimg2++) + (float)*pimg2; 
avg/=8.0; 

*(pgrad++) = (float)*(pimg++) - avg; 

} 

pimg+=2; pimgl +=2; pimg2+=2; pgrad+=2; 

} 

return(O); 
} 



main( argc, argv ) 
Int argc ; 
char *argvQ ; 

{ 

long ij,code_number l go=1,cur_bit,button; 

long temp,incremenUast_middIe,count_thresh = 0; 

long IntJotalJBITS]; 

unsigned char ^pinng.tmp/pimgl *img_out *noise*pnolse; 
char string[80J, outfile[80]; 
FILE *inf; 

float ftemp,current_scale = 1.0,rms,avg; 

float currentjgv = 0.0035, currentjgrad = 0.0027; 

float f_inc,mult; 

float *pgrad,*pgv; 

long gid_img,gid _img1 ,gid_signed,which; 

lf(argc<2 || argc ==3 || argc>4) { 

for(j=0 J<helpiines;j++)fprintf( stderr, "%s M , help[fl) ; 
exit( 1 ) ; 

} 

else if(argc ==2){ 
xdim = 324; 
ydim = 414; 

} 

else { 

xdim = atoi( argv[2] ); 
ydim = atoi( argv[3] ); 

} 

if(DISPLAYJT) disp = calloc(xdim*ydim, sizeof(Co!orindex) ) ; 
else disp = calloc(40000, sizeof(Colorindex) ) ; 

scale_disp = 031100(40000, sizeof(Colorindex) ) ; 

img = calloc(xdim*ydim, sizeof(unsigned char) ) ; 

img_put = calloc(xdim*ydim, sizeof(unsigned char) ) ; 

noise * calloc(xdim*ydim, sizeof(unsigned char) ) ; 

gv = calloc(xdim*ydim, sizeof(float) ) ; 

grad = calloc(xdim*ydim, sizeof(float) ) ; 

if( Idisp || !scate_disp || limg || !img_put || lnoise || !gv || !grad){ 
fprintf( stderr, M sign_it : can not allocate space\n" ) ; 
exit( 1 ) ; 

} 

r read in binary data into array */ 

inf=fopen(argv[1),Y'); 

if(!inf) { 



fprintf(stderr, w sign it: can't open %s\n'\argv[1]); 
exit(1); 

} 

fread(img,sizeof(unsigned charj.xditrfydimjnf); 
fclose(inf); 

rflipit*/ 
pimg = img; 

pimgl = &img[xdim*(ydim-1)]: 
for(i=0;i<(ydim/2);i++){ 

forQ=0;j<xdim-j++){ 

tmp = *pimg; 

*(pimg+.+) = *pimg1; 

*(pimg1++) = tmp; 

} 

pimgl -= (2*xdim); 

} 

r open and read in master noise frame */ 
inf = fopenfmaster noiseYO: 
K(linf) { 

fprintf(stderr f H sign_it: can't open yosXnVmasterjioise*); 
exit(1); 

} 

if((xdim*ydim) > NOISE_MAX ){ 

fprintf(stderr, w sign_it: image size larger than master key M ); 
exit(1); 

} 

fread(noise,sizeof(unsigned char),xdim*ydim,inf); 
fclose(inf); 

r calculate avg then rms noise of masterjioise image */ 
pnoise = noise; 
avg = 0.0; 

for(i=0;i<(ydim*xdim);i++){ 

ftemp = (float)*(pnoise++); 
avg += ftemp; 

) 

avg /= (float)(xdim*ydim); 

printff'AVG noise value of master noise image = %f \n\n*\avg); 
pnoise = noise; 
rms = 0.0; 

for(i=0;i<(ydim*xdim);i++){ 

ftemp = (float)*(pnoise++) - avg ; 
ftemp *= ftemp; 
rms += ftemp; 




} 

rms /= (xdim*ydim); 

rms = (float)sqrt((double)rms); 

printff'RMS noise value of master noise image = %f \n\n",rms); 

r calculate gv */ 

calculate jv(img,gv f xdim l ydim); 

/• calculate grad */ 
calculate^gradfimg.grad.xdlm.ydim); 

r user input */ 
code_number = -1 ; 

white(code_number < 0 || code_number > MAX1){ 

printffType in a number between 0 and %ld: ",MAX1); 
scanf("%d M ,&code number); 

} 

fi!l_bits(code_number); 

if(DISPLAYJT){ 

foreground(); 

prefsize(xdim,ydim); 

gidjmg = winopenfOriginal?"); 

gflush(); 

pdisp = disp; 
pimg = img; 

for(i=0;i<(ydim*xdim);i++){ 

*(pdisp++) = (Colorindex)*(pimg++); 

} 

rectwrite(0,0,xdim-1 ,ydim-1 ,disp); 

prefsize(xdim,ydim); 

gidjmgl = winopenfOrignal?"); 

rectwriteCO.O.xdim-l ,ydim-1 ,disp); 

prefsize(xdim,ydim); 

gid_signed = winopen( ,, Orignal? M ); 

load_scale_disp(); 
prefsize(400 f 50); 

gid_gv = winopen("root GV rms noise"); 
prefsize(400,50); 

gid_grad = winopen("GRAD rms noise"); 

} 
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r main visual feedback loop V 
- last_middle = 0; 
while(go){ 

pdisp = disp; 

pimg = Img; 

pimgl = img_out; 

pnoise = noise; 

pgv = gv; 

pgrad = grad; 

for(i=0;i<(ydim*xdim);i++){ 
cur_bit = i%BITS; 

mult = currentjgv * *(pgv++) + current jgrad * *(pgrad++); 
if( (float)*pnoise > avg ){ 

fjnc = ((float)*(pnoise++) - avg ); 

increment = (int)(f Jnc * mult); 

} 

else { 

fjnc = (avg - (float)*(pnoise++) ); - 
increment = (int)(fjnc * mult); 
increment = -increment; 

} 

if( fabs( (double) *(pgrad - 1)) > GRADJTHRESHOLD ){ 
increment = 0; 

} 

if(bit[cur_bit]){ 

temp = (int)*(pimg) + increment; 

} 

e!se{ 

temp = (int)*(pimg) - increment; 

} 

if( *(pimg++) ~ 0 )tmp=0; 
else if(temp < 0)tmp=0; 
else if(temp > 255)tmp=255; 
else tmp= (unsigned char)temp; 
if(D!SPLAYJT) *(pdisp++) = (Colorindex)tmp; 
*(pimg1++) = tmp; 

} 

if(DISPLAYJT){ 

winset(gid_signed); 
rectwritefO.O.xdim-l ,ydim-1 .disp); 
pdisp = disp; 
pimg = img; 

for(i=0;i<(ydim*xdim);i++){ 

*(pdisp++) = (Colorindex)*(pimg++); 



winset(gidjmg); 

rectwriteCO.O.xdim-l ,ydim-1 ,disp); 
winset(gid_img1); 
rectwrite(dro,xdim-1 ,ydim-1 ,disp); 

draw_scale(rms*current_gv t gid jgv); 
draw_scale(rms*currentjgrad l gidjgrad); 



button=0; 
while(!button){ 

current_scale=1 .0; 
if( getbutton(LEFTMOUSE) ){ 
current_scale 1.15; 
button = 1; 
last_middle = 0; 

} 

if( getbutton(RIGHTMOUSE) ){ 
current_scale *= 0.85; 
button = 1; 
last middle = 0; 

} 

if( getbutton(MIDDLEMOUSE) ){ 

while( getbutton(MIDDLEMOUSE) ); 
if( last_middle ){ 
button =1; 
go=0; 

} 

else{ 

if(which == gid jgv)which=gidjgrad; 
else which=gidjgv; 
Iast_middle = 1 ; 

} 

} 

} ! 
if(go){ I 

if( which " gid jv )current_gv *= current_scale; 

else current_grad *= current_scale; 

} 

} 

r re-calculate and write out master_key */ 
pnoise = noise; 

pimg = img; r use img as a temporary buffer 7 
pgv = gv; 
pgrad = grad; 
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counMhresh = 0; 
for(i=0;i<(ydim*xdim);i++){ 

mult = current jgv * *(pgv++) + currentjgrad * *(pgrad++); 
if( (float)*pnoise > avg ){ 

fjnc « ((float)*(pnoise++) - avg ); 
increment = (int)(fjnc * mult); 

} 

else{ 

Hnc = (avg - (f!oat)*(pnolse++) ); 
increment = (int)(fjnc * mult); 
increment = -increment; 

if( fabs( (double) *(pgrad - 1)) > GRADJTHRESHOLD ){ 
count_thresh++; 
increment = 0; 

} 

*(pimg++) = (unsigned char)128 + (unsigned char)increment; 



r find mean value of the integer based noise, 

for each bit plane, then remove it, and store to dnoise */ 
for(i=0;i<BITS;i++)inUotal[i] = 0; 
pimg = img; 

for(i=0;i<(ydim*xdim);i++){ 

intJotal[i%BITS] += (long) (*(pimg++) - 128 ); 

} 

rms = 0.0; 
for(j=0;j<B!TS;j++){ 

if(intjotai0] > 0){ 

for(i=0;i<inMotalb*];i++){ 

temp = (long)randO; 
temp BITS; 
temp += j; 

if( (temp <(xdim * ydim)) && (img[temp] > 128) 

)img[temp]-; 

else i-; 

} 

} 

else if(intJotalQ] < 0){ 

intjotalfl] = - intjotalfl]; 
for(i=0;i<int_total[j];i++){ 

temp = (long)rand(); 

temp *= BITS; 

temp += j; 

if( temp <(xdim * ydim) && (img[temp] < 128) 

)img[temp]++; 



else i-; 

} 

} 

} 

rms = 0.0; 
pimg = img; 

for(i=0;i<(ydim*xdim);i++){ 

temp * (long) *(pimg++) - 128; 
rms += (float)(temp*temp); 

} 

rms /= (xdim * ydim); 

rms - (float)sqrt((double)rms); 

printfOnFINAL RMS noise value of master key = %f \n\n*\rms); 
printffcount threshold hits = %d \n\n n ,counMhresh); 
sprintf(outfile,"master_key"); 
inf = fopen(outfi!e ( w w M ); 
if(!inf) { 

fprintf(stderr,"sign It: can't open %s\n",outfile); 
exit(1); 

} 

fwrite(img,sizeof(unsigned char) f xdim*ydim,inf); 
fcloseflnf); 

T now re-sign output image because of slight changes to master key */ 



r flip output image */ 
pimg = img_out; 

pimgl = &img_out[xdim*(ydim-1)]; 
for(i=0;i<(ydim/2);i++){ 

for(j=0;j<xdim;j++){ 

tmp = *pimg; 

*(pimg++) = *pimg1; 

*(pimg1++) = tmp; 

> 

pimgl -= (2*xdim); 

} 

r write out signed image */ 
sprintf(outfile ( ,, %s.signed" l argv[1)); 
inf = fopen(outfi!e,"w"); 
if(!inf) { 

fprintf(stderr, M sign_it: can't open %s\n ,, f outfiIe); 
exit(1); 

} 

fwrite(img_out t sizeof(unsigned char),xdim*ydim,inf); 
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fclose(inf); 



r free and clean up */ 
if(DISPLAYJT) gexitO; 
free(img); 
free(img_out); 
free(disp); 
free(sca!e_disp); 
free(noise); 
free(gv); 
free(grad); 



r this program takes in a raw 8 bit image which is in general a sub-section of 
another original image. The orginal image will have already been signed. This 
program needs to know what offest and extent that the sub-section has realtive to 
the dimensions of the original. It also needs to have the master_noise file present 
keyed to the original extent of the image also */ 

r this one now uses the dot product approach to remove the innate bias from the 
original */ 



#include "pinecone.h" 
#include "../.VsgMooIs/image.h" 

#define BITS 16 
#define VALUE 32768 
#define INCREMENT 0.5 

int bit[BITS]; 

long dim,bfts,fftdim,fft_size; 
float *ar,*ai; 

double wr[1 0000],wi[1 0000]; 

int helplines = 7; 
char *heIpQ = 
{ 

"usage: recognize filename originaljfile orig_xdim orig_ydim x_offset y_offset 

patch_xdim patch_jydim\n", 

"\n Program requires that the file 'master-key' exists in current directory", 
"\n filename is image in 8-bit raw form which is in general sub-section of 

original image", 

"\n original_file is original image in 8-bit raw form M , 

"\n orig_xdim and orig_ydim is the pixel extent of the original image", 

"\n x_offset and y_offset is the offset in pixels of the upper left edge of the 

sub-section rel. to original", 

"\n patch_xdim and patch_ydim is the pixel extent of the image patch being 

recognized" 

}; 



int general_purpose_filter(float 'image, long xdim, long ydim, long which){ 
long ij; 

float *par,*pai,*pan,*pan_im,*pimage; 
float k_va!ue,zero_point t iijj,ftemp; 
float std_dev=0.0; 
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memset(ar,0, fft_size * sizeof(f)oat) ); 
memset^ai.O, fft_size * sizeof(float) ); 

r copy image buffer into ar */ 
par = ar; 
pimage = image; 
for(i=0;i<ydim;i++){ 

for(j=0;j<xdinrj++){ 

*(par++) = *(pimage++); 

} 

par += (fftdim - xdim); 

} 

r real 2d fft, in place 7 
reaifft2d(ar,ai,bits t 0,wr ( wi); 

r if which ~ 0 or 1, find the standard deviation in frequency space 7 
if(which == 0 || which == 1 ){ 
par = ar; pai = ai; 

pan = &ar[ffidim-1];panjm=&ai[ffidim-1]; 
std_dev = 0.0; 
for(i=0;i<(1 +fftdim/2);i++){ 

ii = (f!oat)i / (f1oat)(fftdim/2); 

ii *= ii; 

forO=0;j<(1 +fftdim/2);j++){ 

jj = (float)j / (float)(fftdim/2); 
jj = (f]oat)sqrt((double)(ii+ij*ii)); 

ftemp = *par • *par + *pai * *pai; 
std_dev += jj*ftemp; 
par++; pai++; 
ifG!=0&&j!=(fftdim^)){ 

ftemp = *pan * *pan + *panjm * *panjm; 

std_dev += jj*ftemp; 

pan-; panjm-; 

} 

} 

par += (fftdim/2 -1); pai += (fftdim/2 -1); 

pan += (fftdim +fftdim/2 -1); panjm += (ffidim+fftdim/2 -1); 

printf("\n Std deviation found = %e H .std_dev); 



r now filter it in the case of which == 1 7 



if(which == 3){ 

k_value =1.0/ powf(zero_point,3.0); 
tf(zero_point > 0.99 ) zero_point = 0.99; 
par - ar;pai=ai; 

pan = &ar[fftdim-1];pan im=&ai[fftdim-1]; 
for(i=0;i<(1 +ffidim/2);H-+){ 

ii = (float)i / (float)(fftdim/2); 

ii *= ii; 

forG=0;j<(1 +fftdim/2);j++){ 

jj = (float)j / (float)(fftdim/2); 
jj = (float)sqrt((double)(ii+lj*jj)); 

if( jj < zero_point ){ 

jj = 1.0 - (k_value*inj*fl); 

*(par++)*=jj; 

•(pai++)*=ij; 

if(j!=0 &&j!=(fftdim/2)){ 

*(pan-)*=jj; 

*(pan_im-) *= jj; 

} 

> 

else { 

*(par++) = 0.0; 
*(pai++) = 0.0; 
if(j!=0 &&j!=(fttdinV2)){ 

*(pan-) = 0.0; 

*(pan_im-) = 0.0; 

} 

} 

> 

par += (fftdim/2 -1 ); pai += (fftdim/2 -1 ); 

pan += (fftdim +fftdim/2 -1); pan_im += (ffldim+fftdimA2 -1); 

} . 

} 

| 



/* and in the case of which == 2. then use the same filter used in the case 1 */ 



/* in the case of 1 and 2. inverse fft and copy back into the original pointer array */ 



# 



} 
} 

S: return(l): 

r% ' 

*9 



U 
m 



tf( which ~3 || which ==2){ 

realfft2d( ar, ai, bits, 1 , wr, wi) ; 

par = ar; 

pimage = image; 

for(i=0;i<ydim;i++){ 

for(j=OJ<xdim;j++){ 

*(pimage++) = *(par++); 

} 

par += (fftdim - xdim); 



i{| long i; 

3 * float *pimg1 ,*pimg2; 

double result = 0.0; 



double dot_product(float *img1, float *img2,long size){ 



pimgl = imgl; 
pimg2 = img2; 

for(i=0;i<size;i++) result += ( (double) *(pimg1++) * (double) *(pimg2++) ); 

return(result); 
} 

double dot_squared_product(float *img1 .float *img2Jong size){ 
long I; 

float *pimg1 f *pimg2; 
double result = 0.0; 
double dtemp; 

pimgl = imgl; 
pimg2 = img2; 

for(i=0;i<size;i++){ 

dtemp = (double)*(pimg1++); 
if(dtemp < 0.0){ 

dtemp*=dtemp; 
dtemp = -dtemp; 

} 

else dtemp *= dtemp; 
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result += ( dtemp * (double) *(pimg2++) ); 

} 



return(result); 
} 

double dot_root_product(float *img1 .float *img2,tong size){ 
long i; 

float *pimg 1 ,*pimg2; 
double result = 0.0; 
double dtemp; 

pimgl = imgl; 
pimg2 = Img2; 

for(i=0;i<size;i++){ 

dtemp = (double)*(pimg1++); 
rf(dtemp <0){ 

dtemp = sqrt(-dtemp); 
dtemp = -dtemp; 

} 

else dtemp = sqrt(dtemp); 

result += ( dtemp * (double)*(pimg2++) ); 

} 

return(result); 
} 

double dot_rootroot_product(float *img1, float *img2 f long size){ 
long i; 

float *pimg 1 ,*pimg2; 
double result = 0.0; 
double dtemp; 

pimgl = Img1; 
pimg2 = img2; 

for(i=0;i<size;i++){ 

dtemp = (double)*(pimg1++); 
if(dtemp <0){ 

dtemp = sqrt(-dtemp); 
dtemp = sqrt(dtemp); 
dtemp = -dtemp; 

} 

else { 

dtemp = sqrt(dtemp); 
dtemp = sqrt(dtemp); 

} 




1> 



result += ( dtemp * (double)*(pimg2++) ); 

return(result); 
} 



main( argc, argv ) 
int argc; 
char *argvQ ; 
{ 

long ij,k t xdim i ydim,temp l patch - xdim < patch^dim t xoffset i yoffset 1 lndex; 

long value = VALUE,code_number = O.go; 

unsigned char tmp/noise,*pnoise,*original,*poriginal *p1; 

unsigned char *img f *pimg,*pimg1 .strnp; 

char string[80], outfile[80]; 

FlLE*inf; 

float ^fnoise/pfnoise^otal^fimg^pfimg/foriginal/pforiginal; 
float *fnoise_tmp*pfnoise_tmp; 

double base^dot.dotJmg^mag.originaLmag.noise^mag.rms; 

double square_dot,root_dot,rootroot_dot; 

double base_square_dot,base_root_dot,dtemp; 

double imgjTiagj-ooUmg_mag_square ( img_mag rootroot; 

double dot1,root_dot1,root_dot2; 

if(argc != 9) { 

forO=0;j<helplines;j++)fprintf( stderr, ,,0 /oS M , helpm) ; 
exit( 1 ) ; 

} 

xdim = atoi( argv[3] ); 
ydim = atoi( argv[4] ); 
xoffset = atoi( argv[5] ); 
yoffset = atoi( argv[6] ); 
patch_xdim = atoi( argv[7] ); 
patch_ydim = atoi( argv[8] ); 

P find the next power of two equal to or higher than patch */ 
if(patch_xdim > patch_ydim)dim = patch_xdim; 
else dim = patch_ydim; 
fftdim = 1; go = 1; bits = 0; 
while( go ){ 

if( dim > fftdim ){ 

fftdim*=2; 

bits++; 

) 

else go = 0; 

} 



rf(bits>12){ 

fprintf( stderr, -recognize : sorry, this particular program only accepts 
4K images and less\n" ) ; 

exit( 1 ) ; 

} 

ffi_size = fftdim*fftdim; 

img = calloc(patch_xdim*patch_ydim, sizeof(unsigned char) ) ; 
fimg = calloc(patch_xdim*patch_ydim, sizeof(float) ) ; 
original » callocfxdim^ydim, sizeof(uns1gned char) ) ; 
foriginal = ca!loc(patch_xdim*patchjrdim, sizeof(float) ) ; 
noise = calloc(xdim*ydim, sizeof(uns!gned char) ) ; 
fnoise = ca!loc(patch_xdim*patch_ydim, sizeof(float) ) ; 
fnoisejmp = calloc(patch_xdim*patch_ydim, sizeof(float) ) ; 
ar = ca!loc(fft_size, sizeof(float) ) ; 
ai = calIoc(ffi_size, sizeof(float) ) ; 

if( limg || Inoise || loriginal || Ifnoise || !fnoise_tmp || Ifimg || Iforiginal || lar || 

!ai){ 

fprintf( stderr, w sign_it : can not allocate space\n" ) ; 
exit( 1 ) ; 

} 



r read in Suspected Patch into img */ 
inf = fopen(argv[1] f V); 
if(!inf) { 

fprintf(stderr f M signjt: can't open %s\n , \argv[1l); 
exit(1); 

} 

fread( img,sizeof(unsigned char),patch_xdim*patch_ydim,inf); 
fclose(inf); 

/•flip it*/ 
pimg = img; 

p1 = &img[patch_xdim*(patch_ydim-1)]; 
for(i=0;i<(patch_ydim/2);i++){ 

forfl=0;j<patch_xdim;j++){ 

tmp = *pimg; 

*(pimg++) = *p1; 

*(p1++) = tmp; 

} 

p1 -= (2*patch_xdim); 



r copy it into float array, removing mean value */ 
pimg 5= img; 



total = 0.0; 

for(i=0;i<(patch^dim*patch.xdim);i++)total += (float) *(pimg++V 
total /= (f!oat)(patch^xdim*patch _ydim); IP"***)- 
pimg = img; 
pfimg = fimg; 

^T^^'*™-**^ *<P fim 9 ++ > - *(P-9-) - total; 
img_mag = 0.0; 

for(i=0;i<(patch_ydim*patch_xdim);i++K 
dtemp = (double)*(pfimg++); 
dtemp *= dtemp; 
img_mag += dtemp; 

img_mag = sqrt(img_mag); 
printfC , img_mag = %|f ".img_mag); 

filter furTcSon *?* h freqUenCy s P ace usi "9 *• general purpose 

r 



*/ 



generaLpurpose_filter(fimg. P atch_xdim.patch_ydim.O); 

r read in Original Unsigned Image */ 

inf=fopen(argv[2],V); 
if(linf) { 

exit(?) Stderr ' Si9nJt C * rft ° pen %sV, "' ar 9v[2J); 
Jctose^ char ). x dimVdim.inf); 

/* create original patch 7 
index = yoffset * xdim + xoffset; 
poriginal = original; 
p1 = &originalpndex]; 
fora=0;j<patch_ydim;j++){ 

for(k=0;k<patch_xdim;k++){ 

*(poriginal++) = *(p1++); 

^ p1 +■ (xdim - patch_xdim); 

r flip it v 

poriginal = original; 

w i? ri / 9ina, fP atch ^dim*(patch^ydim-1 )J; 
for(i=0;i<(patch_ydirTV2);i+ + ){ 
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for(j ss O;j<patch_xdimj++){ 

tmp = *poriginal; 

*(poriginal++) = *p1; 
^ # (p1++) = tmp; 

^ p1 -= (2*patch_xdim); 

r copy It into float array, removing mean 7 
porigina! = original; 
total = 0.0; 

for(i=0;i<(patch - ydim*patch_xdim);i++)total += (float) *(poriginal++V 
total /= (f!oat)(patch_xdim*patch_ydim); 
poriginal = original; 
pforiginal = foriginal; 

for(i=0;i<(patch^dim*patch_xdim);i++) *(pforigina!++) = (float) *(poriginal++) 

total; 

r now find its std deviation in frequency space, then filter K to roughly match 
the suspect patch 7 
/• 

generaLpurpose_filter(foriginaI,patch_xdim,patch_ydim,1 ); 

pforiginal = foriginal; 
original_mag = 0.0; 

for(i=0;i<(patch_ydim*patch_xdim);i++){ 
dtemp = (double)*(pforigina!++); 
dtemp dtemp; 
original_mag += dtemp; 

original_mag = sqrt(originaLmag); 
printf("originaLmag = %ff \n M ,original_mag); 



r now remove the original image from img by using the dot product or 
whatever 7 

dot = dot_product(forigina!.fimg,patch xdirn*patch_ydim); 
pnntff initial dot = %le H ,dot/(originaLmag*img_mag)); 

dot /= (original_mag * originaLmag); 
printffdot found = %le ",dot); 

pforiginal = foriginal; 

pfimg = fimg; 

img_mag = 0.0; 

for(i=0;i<(patch_ydim*patch_xdim);i++){ 
*pfimg -= (dot * *(pforiginal++)); 



r now re-calculate the magnitudes 7 

if( *pfimg < 0.0 ) dtemp = sqrt( -(doub!e)*pfimg ); 

else dtemp ■ sqrt( (double)*pfimg ); 

img_rnag_rootroot += dtemp; 

dtemp = (doub!enpfimg++); 

if(dtemp < 0.0) img_mag_root -= dtemp; 

else img_mag_root += dtemp; . /■^■■a 
dtemp dtemp; ;r 
img_mag += dtemp; 
dtemp *= dtemp; 
img_mag square += dtemp; 

} . . 

img_mag = sqrt(img_mag); 
img_mag_rootroot = sqrt(img_mag_rootroot); 
img_mag_root = sqrt(img_mag_root); 
img_mag_square = sqrt(img_mag_square); . 
dot = dot_product(foriginaI f fimg,patch_xdim*patch_ydim); 
printff new dot = %le ^dot/(original_mag*img_mag)); 

■ • •••• 'iv^V"' * 

r open and read in master noise frame 7 
inf = fopenCmaster_key",V); 
rf(!inf) { 

fprintf(stderr,"sign_it: can't open %s\nYmaster key"); 
exit(1); 

fread(noise,sizeof(unsigned char) ( xdim*ydim ( inf); ' 
fclose(inf); 

yoffset = ydim - (yoffset + patch_ydim); r needed for flipping 7 



r main recoginition loop 7 

rms = 0.0; 
for(F0;i<BITS;i++){ 

/* create current bit plane noise patch, fnoise 7 
r THE MEAN OF THESE MASTER NOISE PLANES SHOULD BE 
ZERO!!! FROM THE SIGNING PROCESS 7 
index = yoffset * xdim + xoffset; 
pnoise = &noise[index]; 
pfnoise = fnoise; 
forfl=0;j<patch_ydinrj++){ 
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for(k=0;k<patch_xdim;k++){ 
if( (index%BITS) == i ){ 

•pfnoise = (float) (*pnoise - 128); 

rms += (doub!e)*pfnoise * (doub!e)*pfnoise; 

} 

else *pfnoise = 0.0; 

pfnoise++; "\ 

pnoise++; 

index++; 

} 

pnoise += (xdim - patch_xdim); 
index += (xdim - patch xdim); 

} 

r now filter that patch using the same gaussian used on the original V 

r 

general_purpose_filter(fnoise,patch xdim,patch_ydim,2); 

V 

pfnoise = fnoise; 
noise_mag = 0.0; 

forO=0;j<(patch_ydim*patch_xdim)J++){ 

noise_mag += (double) *(pfnoise) * (double) *(pfnolse); 
pfnoise++; 

} 

noise_mag = sqrt(noise_mag); 

r find dot product between suspect patch and filtered tmp noise */ 
dot = dot^roduc^fnoise.fimg.patch^xdim^patch^ydim); 
dot /= (noise_mag * img_mag); 

printf("dot=%le ".dot); 

r new expiriment: re-weight the noise then do the root dot V 
pfimg = fimg; 
pfnoise = fnoise; 
noise_mag = 0.0; 
forG=0;j<patch_ydim;j++){ 

for(k=0;k<patch_xdim;k++){ 

dtemp = (double)*pfimg; 
if( fabs( dtemp ) < 1.0); 
else { 

if( dtemp >= 0.0){ 

•pfnoise /= (floai)sqrt( dtemp ); 

} 

else{ 

. *pfnoise /= (fioat)sqrt( -dtemp ); 
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} 

} 

noisejnag (double) *(pfnoise) * (double) "(pfnoise); 
pfnoise++; pfimg++; 

} 

} 

noise_mag = sqrt(noise_mag); 

root_dot = dot_rootjrodud(fimg.fholse,patch_xdim*patch_ydim); 
root_dot /= (noise_mag * img_mag_root); 
printf( M new root = %le Voot_dot); 

pfimg = fimg; 
pfnoise - fnoise; 
pfnoise_tmp = fnoise_tmp; 
noise_mag = 0.0; 
forQ=0;j<patch_ydim;j++){ 

for(k=0;k<patch_xdim;k++){ 

dtemp = (double)*pfimg; 

if( fabs( dtemp ) < 1 .0 )*pfnoise_tmp = *pfnoise; 

else { 

if( dtemp >= 0.0){ 

*pfnoise_tmp = *pfnoise / (float)sqrt( 

dtemp ); 

} 

else { 

•pfnoise Jmp = *pfnoise / (float)sqrt( - 

dtemp ); 

} 

} 

noise_mag += (double) *(pfnoise_tmp) * (double) 

*(pfnoise_tmp); 

pfnoise++; pfimg++;pfnoise tmp++; 

} 

} 

noise_mag » sqrt(noise_mag); 

root_dot = doLrootjDroductffimg.fnoise.tmp.patch^xdim^patch^ydim); 
root_dot /= (noisejnag * img_mag_root); 
printffnew rootroot = %le Voot_dot); 

P add assumption */ 

r 

pfimg = fimg; 
pfnoise = fnoise; 
pfnoise_tmp = fnoisejmp; 
noise_mag = 0.0; 
forfl=0;j<patch_ydim;j++){ 
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INCREMENT; 



dtemp ); 



dtemp ); 



*(pfnoise_tmp); 



for(k=0;k<patch_xdim;k++){ 

dtemp =7doubIe)*pfimg + (double)*pfnoise * 

lf( fabs( dtemp ) < 1.0 )*pfnoiseJmp = 'pfnoise; 
else{ 

lf( dtemp >= 0.0){ 

*pfrioise_tmp = *pfnoise / (float)sqrt( 



} 

else { 



*pfnoise_tmp = *pfnoise / (float)sqrt( • 



} 



.} 



} 

noise_mag += (double) *(pfnoise Jmp) * (double) 
pfnoise++; pfimg++;pfnoise_tmp++; 



} 

noise_mag = sqrt(noise_mag); 
rooted ot1 = 

dot^rootjDroductCfimg.fnoise^tmp.patch^xdim^patch^ydim); 

root_dot1 /= (no!se_mag * img_mag_root); 
printffnew root add = %le " ( root_dot1); 

r subtract assumption */ 

r 

pfimg = fimg; 
pfnoise = fnoise; 
pfnoise_tmp = fnoise_tmp; 
noise_mag = 0.0; 
forfl=Oj<patch_ydim;j++){ 

for(k=0;k<patch_xdim;k++){ 

dtemp ■» (double)*pfimg - (double)'pfnoise * 



INCREMENT; 



dtemp ); 



dtemp ); 



if( fabs( dtemp ) < 1.0 )*pfnoise_tmp = *pf no ise; 
else { 

if( dtemp >= 0.0){ 

*pfnoise_tmp = *pfnoise / (float)sqrt( 



} 

else { 



•pfnoise Jmp = 'pfnoise / (float)sqrt( - 
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*(pfnoise_tmp); 

} 

} 

noise_mag = sqrt(noise_mag); 
root_dot2 = 

dot_rootj3rodurt(fimg^noise_tmp l patch_xdim*patch_ydim); 

root_dot2 /= (noise_mag * img_mag_root); 
printffnew root sub = %le \n",robt_dot2); 

printf("\ndiff1 = %le diff2 ■ %!e ~,root_dot1 - root_dot,root_dot2- 

root_dot); 



if( dot > 0.0 ){ 
bit[i] = 1; 

code_number+=value; 

} 

else { 

Mt[i]-0; 

} 

printff bit value %d = %d \n H ,i,bit[i]); 

value/=2; 



} 

printf( M \n\nCode Number Found = %d\n\n",code_number); 
rms/=(doubIe)(xdlm*ydim); 
rms = sqrt(rms); 
printfr\n\nrms = VoleXnW'.rms); 



r free and clean up */ 

free(img); 

free(original); 

free(noise); 

free(fnoise); 

free(fnoisejmp); 

free(fimg); 

free(foriginal); 

free(ar); 

free(ai); 

} 



} 

noise_mag += (double) *(pfnoise_tmp) * (double) 
pfnoise++; pfimg++;pfnoise_tmp++; 
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r this program takes in a raw 8 bit image which is in general a sub-section of 
another original image. The orginal image will have already been signed. This 
program needs to know what offest and extent that the sub-section has realtive to 
the dimensions of the original. It also needs to have the master_no1se file present 
keyed to the original extent of the image also */ 

r this one now uses the dot product approach, removing the innate bias from the 
original */ 



#include "pinecone.h" 
#include "../../sgMools/image.h" 

#def?ne BITS 16 
#define VALUE 32768 
#define INCREMENT 0.5 
#define KCONST 500.0 
#define ZERO_POINT 0.2 

int bit[BITS]; 

long dim.bits.fftdim.fftjsize; 
float *ar,*ai; 

double wr[10000],wi[10000]; 

int helplines = 7; 
char *helpQ = 
{ 

"usage: recognize filename original_file orig_xdim orig_ydim x offset y offset 

patch_xdim patch_ydim\n", 

"\n Program requires that the file 'master-key 1 exists in current directory", 
"\n filename is image in 8-bit raw form which is in general sub-section of 

original image", 

"\n priginal_file is original image in 8-bit raw form ", 

"Vn orig_xdim and orig_ydim is the pixel extent of the original image", 

"\n x_offset and y_offset is the offset in pixels of the uppe'r left edge of the 

sub-section rel. to original", 

"\n patch_xdim and patch_ydim is the pixel extent of the image patch being 
recognized" 



int general_purposeJilter(float 'image, long xdim, long ydim, long which){ 
longij; 

float *par,*pai t *pan,*panjm,*pimage; 
float *par1 ,*pai1 t *pan1 ,*pan _im1 ; 
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Appendix D 



float k_va!ue.ii jj f ftemp,ftemp1 .ratio; 
float std_dev, total; 
static float patch_std_dev; 
static float original_std_dev; 

memset(ar,0, fft_size * sizeof(float) ); 
memset(ai,0, fft_size * sizeof(float) ); 

r copy image buffer into ar */ 
par = ar; 

pimage = image; 
for(i=0;i<ydim;i++){ 

for(j=Oj<xdim*j++){ 

*(par++) = *(pimage++); 

} 

par += (fftdim - xdim); 

} 

r real 2d fft, in place */ 
realffEdfar.ai.bits.O.wr.wi); 

r apply low pass to everything */ 
par = ar;pai=ai; 

pan = &ar[fftdim-1];panjm=&ai[fftdim-l]; 
pari = &ar[fft_size - fftdimj;pai1=&ai[fft_size-fftdim]; 
pan1 = &ar[fft_sfce-1];panjm1=&ai[fft size-1]; 
for(i=0;i<(1 +fftdinrV2);i++){ 

ii = (float)i / (fIoat)(fftdim/2); 

ii *= ii; 

for(j=0;j<(1 +fftdim/2);j++){ 

jj = (float)j / (float)(fftdim/2); 
ij = (float)sqrt( (double)(ii+rii) ); 
if(jj>ZERO_POINT)ij=0.0; 
else{ 

jj = (ZEROJ>OlNT- jj)/ZERO_POINT; 

} 

*par*=jj; 

*pai *=ij; 

par++; pai++; 

ifG!=0 &&j!=(fftdim/2)){ 

*pan *=jj; 

*pan Jm *= jj; 

pan-; panjm-; 

} 

} 



# 



par += (ffidim/2 -1 ); pai += (fftdim/2 -1 ); 
pan += (ffidim +fftdim/2 -1); pan_im += (fftdim+ffidim/2 -1); 

if(i!=0 && H=(fftdim/2)){ 

for(j=0;j<(1 +ffldim/2)j++){ 

jj = (float)j / (float)(fftdfnV2); 
jj = (float)sqrt( (double)(ii+jj*jj) ); 
if(jj > 2ERO_POINT )jj =0.0; 
else { 

} j» = (ZERO.POINT - jj)/ZERO_POINT; 

*(par1++)*=ii; 
*(pai1++)«=jj; 
ifONO &&j!=(fftdinV2)){ 
*(pan1-)*=jj; 
•(pan_im1-.)*=jj; 

H } 

Wi P ar1 (3*fftdim/2 +1); pail -= (3*fftdim/2 +1); 

Pan1^(fftdirtV2+1);pan_im1-=(ffidinV2+1); 



Hi 



> 



SSSKL" ; JSiS 5f Standafd deViati0 " in *»*""» space v 
O par « ar; pai * ai; 

W P 3 " = &ar I«dim-1];paruim=&ai[fftdim-1l; 

std_dev = 0.0; 
total = 0.0; 

for(i=0;i<(l +fftdim/2);i++){ 

» s (float)i/(float)(fftdim/2); 
ii *= ii; . 

for(j=0;j<(1 +fftdim/2);j++){ 

jj = (float)j / (float)(ffidim/2); 
jj = (float)sqrt((double)(ii+jj*|j)); 

ftemp = *par * *par + *pai * *pai; 

std_dev += jj*ftemp; 

total += ftemp; 

par++; pai++; 

ifO!=0 &&jl=(ffidim/2)){ 

ftemp = -pan * -pan + -panjm * *pan im; 

total += ftemp; " 

std_dev +=jj*ftemp; 

pan-; panjm-; 
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) 



} 

) 

par += (fftdim/2 -1); pai (fftdim/2 -1); 

pan += (ffidim +fftdim/2 -1); panjm += (fftdim+ffldim/2 -1); 

std_dev/= total; 

printffVn Std deviation found = %e ",std_dev); 



if(which == 0) patch_std_dev = std_dev; 
if(which == 1) original_std_dev = std_dev; 



r now filter it in the case of which ~ 1 */ 

r and in the case of which ~ 2, then use the same filter used in the case 1 7 
if( (which == 1 || which == 2) && (original_std_dev > 0.95 * patch_std_dev) ){ 
ratio = patch_std_dev / original_std_dev; 
ratio ; 

if(ratio > 1.0)k_value = 0.0; 
e!se{ 

k_value = -KCONST * (1 .0 - ratio) ; 

std_dev /= total; 
par = ar;pai=ai; 

pan = &arlfftdim-1];panjm=&ai[fftdim-1]; 
pari = &ar[ffi_size - ffidim];pai1=&ai[ffi_stee-fftdim]; 
iW Pan1 = &ar[fft_size-1];pan im1=&ai[fft size-1]; 

lip std_dev = 0.0; 

total = 0.0; 
< for(i=0;i<(1+fftdim/2);i++){ 

ii = (float)i / (fioat)(ffidim/2); 
ii *= ii; 

for(j=0;j<(1 +fftdim/2);j++){ 

jj = (float)j / (fioat)(fftdim/2); 
jj - ii+jj*jj; 

ftempl = (float)sqrt((double)lj ): 

jj = fexp(jj*k_value); 

*par*=jj; 

*pai *= jj; 

ftemp = *par * *par + *pai * *pai; 
std_dev += ftempl *ftemp; 
total += ftemp; 
par++; pai++; 
if(j!=0 && j!=(fftdim/2)){ 
*pan*=jj; 



*panjm *» jj; 

ftemp = *pan * *pan + *panjm * *pan_im; 
std_dev += ftemp1*ftemp; 
total += ftemp; 
pan-; panjm-; 

} 

} 

par += (fftdim/2 -1); pai += (fftdim/2 -1); 

pan += (ffidim +fftdim/2 -1); panjm += (fftdirn+fftdim/2 -1); 

if(i!=0 && i!=(fftdim/2)){ 

forG=0;j<(1+ffidim/2)J++){ 

jj = (float)j / (float)(fftdim/2); 
jj = ii+iTJj; 

jj = fexp(jj * k_value); 
*(part++) *=jj; 
*(pai1++)*=ij; 
rfQ!=0&&j!=(fftdinV2)){ 

*(pan1-)*=jj; 

*(panjm1-)*=jj; 

} 

} 

pari -= (3*fftdim/2 +1); pail -= (3*fftdim/2 +1); 
pan1 -= (fftdim/2 +1); pan im1 -= (fftdim/2 +1); 

} 

} 

std_dev/= total; 

if(which == 1) printf("\n Std deviation post-filtering = %e M ,std_dev); 
} 



realfft2d( ar, ai, bits, 1 , wr, wi) ; 
par = ar; 

pimage = image; 
for(i=0;i<ydim;i++){ 

for(j=0;j<xdim;j++){ 

*(pimage++) = *(par++); 

} 

par += (fftdim - xdim); 

} 

return(1); 



• 



} 



double doU>roduct(float *img1, float *img2,long size){ 
long i; 

float *pimg1 *pimg2; 
double result = 0.0; 



pimgl = imgl; 
P pimg2 = img2; 



O for(i=0;i<size;i++) result += ( (double) *(pimg1++) •* (double) *(pimg2++) ); 



n > 

Nl 



return(result); 



double dtemp; 



double dot_squared_product(float *img1, float *img2,long size){ 
long i; 
float *pimg1 t *pimg2; 
double result = 0.0; 



pimgl = imgl; 
pimg2 = img2; 
j;| for(i=0;i<size;i++){ 
m dtem P = (doubie)*(pimg1++); 

if(dtemp < 0.0){ 

dtemp*=dtemp; 
dtemp = -dtemp; 

else dtemp *= dtemp; 
. result += ( dtemp * (double) *(pimg2++) ); 



return(result); 
} 

double dot_root_product(float *img1, float *img2,Iong size){ 
long i; 

float *pimg1.*pimg2; 
double result = 0.0; 
double dtemp; 

pimgl = imgl; 
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pimg2 = img2; 

for(i=0;i<size;i++){ 

dtemp = (double)*(pimg1++); 
if(dtemp <0){ 

dtemp = sqrt(-dtemp); 
dtemp = -dtemp; 

} 

else dtemp = sqrt(dtemp); 

result += ( dtemp * (double)*(pimg2++) ); 

retum(result); 
} 

double dotjrootroot_product(float *img1 .float *img2,long size){ 
long i; 

float *pimgVpimg2; 
double result = 0.0; 
double dtemp; 

pimgl = imgl;* 
pimg2 = img2; 

for(i=0;i<size;i++){ 

dtemp = (double)*(pimg1++); 
lf(dtemp <0){ 

dtemp = sqrt(-dtemp); 
dtemp = sqrt(dtemp); 
dtemp = -dtemp; 

} 

else { 

dtemp = sqrt(dtemp); 
dtemp a sqrt(dtemp); 

} 

result += ( dtemp * (double)*(pimg2++) ); 

return(result); 



main( argc, argv ) 
Int argc ; 
char *argvQ ; 
{ 

long ij.k.xdim.ydim.temp.patch^xdim.patch^dim.xoffset.yoffset.index; 
long value = VALUE,code_number = 0,go; 
unsigned char tmp p *nolse f *pnoise ( *original ( *poriginal f *p1; 
unsigned char *img,*pimg,*pimgi ,stmp; 



m 



charstring[80], outfile[80]; 
FILE *inf; 

float ^noise^pfnoise.total/fimg^pfimg^foriginal/pforiginal; 
float *fnoise_tmp *pfnoise Jmp; 

double base_dot,dot l img_mag 1 originaLmag l noise_mag,rms: 

double square_dot l root_dot,rootroot_dot; 

double base_square_dot,basej-oot_dot,dtemp; 

double »mg_mag_root i lmg_mag_square l img_mag_rootroot; 

double dot1 , root_dot1 t root_dot2; 

if(argc l= 9) { 

forfl=0;j<helpHnes;j++)fprintf( stderr, *%s", helpfl]) ; 
exit( 1 ) ; 

} 

xdim = atoi( argv[3] ); 
ydim = atoi( argv[4] ); 
xoffset = atoi( argv[5] ); 
yoffset = atoi( argv[6] ); 
patch_xdim = atoi( argv[7] ); 
patch_ydim = atoi( argv[8] ); 

r find the next power of two equal to or higher than patch */ 
if(patch_xdim > patch_ydim)dim - patch_xdim; 
else dim = patch_ydim; 
fftdim = 1; go * 1; bits = 0; 
whlle( go ){ 

if( dim > fftdim ){ 

ffidim*=2; 

bits++; 

j 

else go = 0; 
if(bits>12){ 

4K images anX^T' ™* 9 "* ° n * 3CCeptS 

exit( 1 ) ; 

fft_size = fftdim * fftdim; 

img = calloc(patch_xdim*patch_ydim. sizeof(unsigned char) ) ■ 
fimg = calloc(patch_xdim*patch_ydim, sizeof(float) ) ; 
original = calloc(xdim*ydim. sizeof(unsigned char) ) •' 
fongmal = calloc(patch_xdim*patch_ydim, sizeof(float) ) • 
noise = calloc(xdim*ydim, sizeof(unsigned char) ) ■ 
fnoise = calloc(patch_xdim*patch_ydim. sizeof(float) ) ■ 
fnoisejmp = calloc(patch_xdim*patch_ydim, sizeof(float) ) • 
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ar « callocCfft^size, sizeof(float) ) ; 
ai = calloc(fft_size ( sizeof(float) ) ; 

if( limg || Inoise || loriginal || Ifnoise || Ifnoise_tmp || Ifimg || Iforiginal || lar || 

fprintf( stderr, "sign_it : can not allocate spaceVT ) ; 
exit( 1 ) ; 

} 



r read in Suspected Patch into img V 
inf = fopen(argv[1] ( V); 
ifflinf) { 

fprintf(stderr,"sign_it: can't open %s\n M ,argv[1]); 
exit(1); 

} 

fread( img .sizeof (unsigned char),patch_xdim*patch_ydim t inf); 
fclose(inf); 

r flip it*/ 
pimg = img; 

p1 = &img[patch_xdim*(patch_ydim-1)]; 
for(i=0;i<(patch_ydim/2);i++){ 

for(j=0;j<patch_xdimJ++){ 

tmp = *pimg; 

*(pimg++) = *p1; 

*(p1++) = tmp; 

} 

p1 -= (2*patch_xdim); 

} 

r copy it into float array, removing mean value */ 
pimg = img; 
total a 0.0; 

for(i=0;i<(patch^ydim*patch_xdim);i-H-)totaI += (float) *(pimg++); 
total /= (float)(patch_xdim*patch_ydim); 
pimg = img; 
pfimg = fimg; 

for(i=0;i<(patch_ydim*patch_xdim);i++) *(pfimg++) = (float) *(pimg++) - total; 
pfimg = fimg; 
img_mag = 0.0; 

for(i=0;i<(patch_ydim*patch_xdim);i++){ 
dtemp = (double)*(pfimg++); 
dtemp *= dtemp; 
img mag += dtemp; 

> 

img_mag = sqrt(img_mag); 



printf("img_mag = %lf Mmgjnag); 

/* now find its std deviation in frequency space using the general purpose 
filter function */ 

general^purpose^filtertftmg.patch^xdim.patch^ydim.O); 

r read in Original Unsigned Image V 
inf = fopen(argv[2],V); 
if(linf) { 

fprintf(stderr, M sign_it: can't open %s\n",argv[2]); 
exit(1); 

} 

fread(original,sizeof(unsigned char),xdim*ydim,inf); 
fclose(inf); 

f create original patch */ 
index = yoffset * xdim + xoffset; 
poriginal « original; 
p1 = &original[index]; 
for(j=0;j<patch_ydim;j++){ 

for(k=0;k<patch_xdim;k++){ 

*(poriginal++) ■ *(p1++); 

} 

p1 += (xdim - patch_xdim); 

} 

r flip it */ 

poriginal = original; 

p1 = &original[patch_xdim*(patch_ydim-1)]; 
for(l=0;i<(patch_ydim/2);i++){ 

for(j=Oj<patch_xdimy++){ 

tmp = 'poriginal; 

*(poriginal++)s*p1; 

*(p1++) = tmp; 

} 

p1 — (2*patch_xdim); 

} 

1* copy it into float array, removing mean */ 
poriginal = original; 
total = 0.0; 

for(i=0;i<(patch_ydim # patch_xdim);i++)total += (float) *(poriginal++); 
total /= (float)(patch_xdim*patch_ydim); 
poriginal = original; 
pforiginal = foriginal; 



10 





J 



for(i=0;i<(patch^dim*patch_xdim);i++) *(pforiginal++) = (float) *(poriginal++) 

total; 

r now find its std deviation in frequency space, then filter it to roughly match 
the suspect patch */ 

general^purpose^filtertforiginal.patch^xdim^atch^ydimj); 

pforiginal = foriginal; 
originaLmag = 0.0; 

for(i=0;i<(patch_ydim*patch_xdim);i++){ 
dtemp = (double)*(pforiginal++); 
dtemp *= dtemp; 
originaLmag +« dtemp; 

} 

originaLmag = sqrt(original_mag); 
printfforiginaLmag = %lf \n" f originaLmag); 



r now remove the original image from img by using the dot product, or 
whatever */ 

dot = doLproductCforiginal.fimg.patch^xdim^patch^ydim); 
printff initial dot = %le ".dot/toriginaLmag^mg^mag)); 

dot /= (originaLmag * originaLmag); 
printffdot found = %le ",dot); 

pforiginal = foriginal; 

pfimg = fimg; 

img_mag = 0.0; 

for(i~0;i<(patch_ydirirpatch_xdim);i++){ . 



*pfimg -« (dot * *(pforiginal++)); 

r now re-calculate the magnitudes */ 

if( *pfimg < 0.0 ) dtemp ■ sqrt( -(double)*pfimg ); 

else dtemp = sqrt( (double)*pfimg ); 

img_mag_rootroot += dtemp; 

dtemp = (double)*(pfimg++); 

if(dtemp < 0.0) img_mag_root -= dtemp; 

else img_mag_root += dtemp; 

dtemp *= dtemp; 

img_mag += dtemp; 

dtemp *= dtemp; 

img_mag_square += dtemp; 



} 

img_mag = sqrt(img_mag); 
img_mag_rootroot = sqrt(img_mag_rootroot); 
imgjnag_root = sqrt(img_mag^root); 
irhg_mag_square = sqrt(img_mag_square); 



dot = dot_product(foriginal ( fimg,patch_xdim*patch_ydim); 
printff new dot = %le ^dot/(originaLmag*img_mag)); 



r open and read in master noise frame */ 

inf = fopenfmasterJceyYO; 
if(!inf) { 

fprintfCstderr.-signJt: can't open yostaVmasterJcey"); 
exit(1); 

> 

fread(noise,sizeof(unsigned char),xdim*ydimjnf); 
fclose(inf); 

yoffset = ydim - (yoffset + patch_ydim); f needed for flipping */ 



r main recoginition loop */ 

rms = 0.0; 
for(i=0;i<BITS;i++){ 

r create current bit plane noise patch, fnoise */ 
r THE MEAN OF THESE MASTER NOISE PLANES SHOULD BE 
ZERO!!! FROM THE SIGNING PROCESS 7 
index = yoffset * xdim + xoffset; 
pnoise = &noise[index]; 
pfnoise = fnoise; 
for(j=0;j<patch_ydimy-H»){ 

for(k=0;k<patch_xdim;k++){ 
rf((index%BITS)==i){ 

*pfnoise = (float) (*pnoise - 1 28); 

rms += (double)*pfnoise * (doub!e)*pfnoise; * 

else *pfnoise = 0.0; 
pfnoise++; 
pnoise++; 
index++; 

> 

pnoise += (xdim - patch_xdim); 
index += (xdim * patch_xdim); 

} 

I* now filter that patch using the same gaussian used on the original V 

generaLpurpose^filte^fnoise.patch^xdim.patch^ydim^); 
pfnoise = fnoise; 
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noise_mag = 0.0; 

for(j=0;j<(patch _ydim*patch_xdim)J++){ 

noise_mag (double) *(pfnoise) * (double) '(pfnoise); 
pfnoise++; 

} 

noise_mag = sqrt(noise_mag); 

r find dot product between suspect patch and filtered tmp noise 7 
dot = dot_product(fnoise 1 fimg,patch_xdim*patch_ydim); 
dot /= (noise_mag * img_mag); 

printf("dot=%le M ,dot); 

r new expiriment: re-weight the noise then do the root dot 7 
pfimg = fimg; 
pfnoise = fnoise; 
noise_mag = 0.0; 
for(j=0;j<patch_ydim;j++){ 

for(k=0;k<patch_xdim;k++){ 

dtemp = (double)*pfimg; 
if( fabs( dtemp ) <1.0); 
else { 

rf( dtemp >= 0.0){ 

'pfnoise /= (float)sqrt( dtemp ); 

} 

else { 

'pfnoise /= (float)sqrt( -dtemp ); 

} 

} 

noise_mag += (double) '(pfnoise) * (double) ^(pfnoise); 
pfnoise++; pfimg++; 

} 

. } 

noise_mag = sqrt(noise_mag); 

root_dot = dot^root^productffimg.fnoise.patch^xdim^patch^ydim); 
root_dot /= (noise_mag * img_mag_root); 1 
printff new root = %!e Voot_dot); 

pfimg = fimg; 
pfnoise = fnoise; 
pfnoise_tmp = fnoise_tmp; 
noise_mag = 0.0; 
for(j=0;j<patch_ydim;j++){ 

for(k=0;k<patch_xdim;k++){ 

dtemp = (double)*pfimg; 

if( fabs( dtemp ) < 1.0 J'pfnoisejmp = 'pfnoise; 



K \1 



ft 



dtemp ); 



dtemp ); 



*(pfnoise_tmp); 



lf( dtemp >= 0.0){ 

*pfnoise_tmp = *pfnoise / (float)sqrt( 

} 

else { 

*pfnoise_tmp = *pf n o'se / (float)sqrt( • 

} 

} 

noise_mag += (double) *(pfnoise_tmp) * (double) 



pfnoise++; pfimg++;pfnoise_tmp++; 
} 

} 

noise_mag = sqrt(noise_mag); 

root_dot = dot^rootjrodurttfimg.fnoise^tmp.patch^xdim^patch^dim); 
root_dot/= (noise_mag * img_mag_root); 
printffnew rootroot = %le Voot_d"ot); 

/* add assumption */ 

r 

pfimg = fimg; 
pfnoise = fnoise; 
pfnoise_tmp = fnoise_tmp; 
noise_mag = 0.0; 
forO=OJ<patch_ydim;j++){ 

for(k=0;k<patch_xdim;k++){ 

dtemp = (double)*pfimg + (double)*pfnoise * 



INCREMENT; 

dtemp ); 

dtemp ); 

*(pfnoise_tmp); 

} 



if( fabs( dtemp ) < 1.0 )*pfnoise_tmp = *pfno!se; 
else { 

if( dtemp >= 0.0){ 

*pfnoise_tmp = *pfnoise / (float)sqrt( 

) 

else{ 

•pfnoise Jmp = *pfnoise / (float)sqrt( - 

} 

) 

noise_mag += (double) *(pfnoise_tmp) • (double) 
pfnoise++; pfimg++;pfnoise_tmp++; 



• ft. 



noise_mag « sqrt(noise mag); 
root_dot1 = 

doLroot^producttfimg.fnoise.tmp.patch^xdim^patch^dim); 

root_dot1 /= (noise_mag * img mag root); 
printffnew root add ■ %Ie w ,root_dot1); 

r subtract assumption 7 

r 

pfimg = fimg; 
pfnoise = fnoise; 
pfnoisejmp = fnoise_tmp; 
noise_mag = 0,0; 
forO=Oj<patch_ydim;j++){ 

for(k=0;k<patch_xdim;k++){ 

INCREMENT* dt6mP = (doub,e)#pfim 9 " Wouble)*pfnoise * 

if( fabs( dtemp ) < 1.0 )*pfnoise tmp = *pfnoise- 
else { . 
if( dtemp >= 0.0){ 

dtemp ); *pfnoise Jmp = *pfnoise / (ftoat)sqrt( 

} 

else { 

dtemp ); fnoise Jmp = 'pfnoise / (float)sqrt( - 

} 

> 

•(pfnoisejmp); "oise^mag += (double) -(pfnoisejmp) * (double) 

pfnoise++; pfimg++;pfnoise_tmp++; 

• } 

noise_mag = sqrt(noise_mag); 
root_dot2 = 

doLrootprodurtffimg.fnoiseJmp.patch.xdim-patch^dim); 

root_dot2 /= (noise_mag * img mag_joot); 
printffnew root sub = %le \nVoot dot2)- 

root.dot); print ^ dlff1 = %,e diff2 = %te Vootjdotl - rooLdot,roo L dot2- 



7 



if( dot > 0.0 ){ 
bHflJ-1; 

code_number+=value; 



} 

else{ 

bit[i] = 0; 

} 

printff bit value %d * %d \n H f i t bitni); 

value/=2; 



printf("\n\nCode Number Found = %d\n\n*,code_number); 
rms/=(double)(xdim*ydim); 
rms = sqrt(rms); 
printf("\n\nrms = %le\n\n\rms); 



r free and clean up */ 

free(img); 

free(original); 

free(noise); 

free(fnoise); 

free(fnoise_tmp); 

free(fimg); 

free(foriginal); 

free(ar); 

free(ai); 
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